International Journal of Advanced Engineering, Management and Science (IJAEMS) 
https://dx. doi. ora/10.22161/ijaems.4.1 .7 


[Vol-4, Issue-1, Jan- 2018] 
ISSN: 2454-1311 


Simulation modeling of the sensitivity analysis of 
differential effects of the intrinsic growth rate of a 
fish population: its implication for the selection of 

a local minimum 

Nwachukwu Eucharia C. 1 , Ekaka-a Enu-Obari N. 2 , Atsu Jeremiah U . 3 


'Department of Mathematics and Statistics, University of Port Harcourt, Port Harcourt, Nigeria 
department of Mathematics, Rivers State University, Nkpolu.Port Harcourt, Nigeria 
"Department of Mathematics/Statistics, Cross River University of Technology, Calabar, Nigeria. 


Abstract — The vulnerability of the differential effects of the 
intrinsic growth rates of the fish population on the 
uncertainty analysis can only be controlled by using the 
mathematical technique of a sensitivity analysis that is 
called a local minimum selection method based on a Matlab 
numerical scheme of ordinary differential equations of 
order 45 (ODE 45). The quantification of the p-nonns 
sensitivity analysis depends on the application of the 1- 
nonn, 2-norm, 3-norm, 4-norm, 5-norm, 6-norm and 
infinity-norm. In the context of this study, the best-fit 
intrinsic growth rate of fish population with a small error 
has occurred when its value is 0.303 which minimizes the 
bigger sensitivity values previously obtained irrespective of 
the p-norm sensitivity values. The novel results which we 
have obtained have not been seen elsewhere. These results 
are fully presented and discussed in this study. 

Keywords — Uncertainty analysis, differential effects, p- 
norms sensitivity analysis, intrinsic growth rate, local 
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I. INTRODUCTION 

Following Ekaka-a et al (2012), sensitivity analysis is a 
measure of the resultant effect due to a variation of a model 
parameter value on the model solution trajectories. It a 
mathematical tool to enhance model validation and 
prediction. On the other hand, on-going debate among many 
researchers including Palumbi (1999), Sumali (2002), 
Pitchford (2007) and Bohnsack (1993) has shown that 
marine reserve is not only economically important but can 
also serve as a tool for equitable management of 
biodiversity especially in the context of fisheries.Kar and 
Charkraborty (2009), in their work titled marine reserves 
and its consequences as a fisheries management tool 


described a prey-predator type fishing model with prey 
dispersed in a 2-patch environment, one of which is called a 
free fishing zone and another, a protected zone. Their main 
method of investigation uses the simulation process. One 
key contribution from their work states that prey-predator 
dichotomy do not matter when implementation of a reserve 
is considered. Their second result shows that reserves will 
be most effective when coupled with fishing effort controls 
in adjacent fisheries. Despite the fact that marine reserves 
and its consequences can be effectively utilized as a 
fisheries management tool, it is still an open research 
problem that these authors did not consider the technique of 
sensitivity analysis which is vital numerical incentive in a 
decision process that can lead to an effective fisheries 
management. 

It remains an open problem to study the differential effects 
of varying the intrinsic growth rate of the fish population on 
the uncertainty analysis using a one-at-a-time sensitivity 
analysis Hamby 1995. It is against this background that we 
propose to use ODE45 RungeKutta numerical scheme with 
initial condition 2 4 2 over a period of fifty (50) weeks to 
study the differential effects of the intrinsic growth rate of 
fish population on the sensitivity analysis that is indexed by 
seven classifications of the sensitivity analysis, namely: 1- 
norm error analysis, 2-norm error analysis, 3-norm error 
analysis, 4-norm error analysis, 5-norm error analysis, 6- 
norm error analysis and infinity-norm error analysis. 

II. MATHEMATICAL FORMULATIONS 

The mathematical model for this research work is based on 
the published article by Kar and Charkraborty (2009) given 
by the first order non-linear ordinary differential equations 


www.iiaems.com 


Page | 29 





International Journal of Advanced Engineering, Management and Science (IJAEMS) [Vol-4, Issue-1, Jan- 2018] 

https://dx. doi. ora/10.22161/ijaems.4.1.7 ISSN: 2454-1311 


which describes the prey-predator fisheries management 
model. Hence, we modify the model as follows: 


^x(t) = ?ix(t) - /?i*(t) 2 - /Hx(t) y(t) - ~x(t) + 
~y(0 ~ mx(t)z{t) (3.1) 

^y(t) = r 2 y(t) - P 2 x(t)y{t ) - /r 2 y(t) 2 + f *(t) - 
~y(0 - ny(t)z(t) - qEy(t ) (3.2) 


^z(t) = r 3 z(t) - 


7~ 3 rzft) 2 

*(t)+y(f) 


(3.3) 


x(0) = 2,y(0) = 4, z(0) = 2 


(3.4) 


In this model, x(f) represents the fish stock in the reserved 
area at time f, y(t) is the fish stock in the unreserved area at 
time t, z(t) is the biomass density of the predator species at 
time t. rys the intrinsic growth rate of fish stock x, r 2 is the 
intrinsic growth rate of fish stocky, r 3 is the intrinsic growth 
rate of the predator species, y is the equilibrium ratio 
between prey biomass and predator biomass, m is maximum 
relative increase in predation to the reserved area or simply 
put, the contribution of z to inhibit the growth of fish stock 
(x) in the reserved area, n is maximum relative increase in 
predation to the unreserved area or simply put, the 
contribution of z to inhibit the growth of fish stock (y) in the 
unreserved area. The predation terms are therefore defined 
as mxz and nyz with respect to fish stocks in the reserved 
and unreserved area respectively. In addition, o is the 
mobility coefficient, a is the size of the reserved area, 1 — 
a is the size of the unreserved area, q is catch ability 
coefficient, E is effort applied for harvesting fish population 
in the unreserved area. 

Further mathematical interpretation can be invoked to 
describe the interaction of the three equations. In equation 
(3.1), — x(t) = rjX(t) shows that without interaction with 

other y and z species, the x population grows unboundedly 
as time increases. The same observation is made for the y 
and z populations. This is mathematically correct but 
practically unrealistic, hence the need for the interaction 
between the species. We can deduce the following 
interpretations term by term: 

^x(t) 2 is the contribution of fish stock in the reserved area 
to inhibits its original exponential growth. 

/r 1 x(t)y(t)is the contribution of fish stock in the 
unreserved area to inhibit the growth of fish population in 
the reserved area. 


-x(t)is the effect of the ratio of the migration rate/mobility 

coefficient to the size of the reserve area to inhibit the 
growth of fish stock in the reserved area 
-^-j-y(t)is the effect of the ratio of the migration 
rate/mobility coefficient to the size of the unreserved area to 
enhance the growth of fish stock in the reserved area 
mx(t)z(t )is the contribution of z to inhibit the growth of 
fish stock on the reserve area 

fi z x(t)y(t)\s the contribution of fish stock in the reserved 
area to inhibit the growth of fish population in the 
unreserved area. 

Ii 2 y(t) 2 \s the contribution of fish stock in the unreserved 
area to inhibits its original exponential growth. 

-x(t)is the effect of the ratio of the migration rate/mobility 

coefficient to the size of the reserve area to inhibit the 
growth of fish stock in the unreserved area 
-^-j-y(t)is the effect of the ratio of the migration 
rate/mobility coefficient to the size of the unreserved area to 
inhibit the growth of fish stock in the unreserved area 
ny(t)z(t) is the contribution of the predator biomass to 
inhibit the growth of fish stock on the unreserved area 
qEy(t) is the contributed effect of the fishing effort and 
Catchability coefficient to inhibit the growth of the fish 
stock in the unreserved area. 

r 3 r z (t) - s equilibrium ratio on the fish stocks in the 

x(t)+y(t) 

reserved and unreserved area due to activities of the 
predator biomass. 

For the purpose of our mathematical analysis, the value for 
model parameter r is the same value for: r 1( r t , /? 1; /? 2 , 
/r 1 ,/r 2 and s — r 3 . In addition, we shall adopt the model 

parameter values as proposed by Kar and Charkraborty 

(2009): r = 0.3, a — 0.2, m = 1, n = 1, q = 

0.0025, s = 0.1, y =0.95 ,E = 95. 

III. METHOD OF ANALYSIS 

The core part of the algorithm which we have utilized to 
calculate the sensitivity of a model parameter is hereby 
described by the following steps which has also been 
implemented in the work of Ekaka-a et al (2013): 

Step I: identify and code the control system of given 
model equations of continuous non-linear first order 
ordinary differential equation in which the model parameter 
is not varied. For the purpose of this analysis, the three 
solution trajectories are denoted by x,y, and z. 

Step II: identify and code a sub-model of the control 
system of given model equations of continuous non-linear 
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first order ordinary differential equation in which the model 
parameter is varied one-at-a-time. In this case, the three 
solution trajectories are denoted by x m ,y m , and z m . 

Step III : code an appropriate Matlab program using ODE 
Runge-Kutta scheme to execute the program in Step I and 
Step II. With the initial conditions and a time range, the 
execution program will produce the solution trajectories for 
the programs in step I and step II. 

On the execution program, specify the difference between 
the solution trajectories of the codes in step I and step II 
as F 1 = x - x m . F 2 — y — y-m and F 3 = z - z m 
Step IV: Use the execution program to calculate the 1- 
norm, 2-norm, 3-norm, 4-norm, 5-norm, 6 -norm and 
infinity-norm for the three solution trajectories of the 
control model equations and similarly for the solution 
trajectories for the difference between the solution 
trajectories. For example, for the x and x m solution 
trajectories which assume precise data points such as Xj and 
x m j, where the subscript j takes on the values of 1, 2, 3, 4, 
...,n, the 1 -norm for the x solution trajectory is defined as 
the sum of the absolute values of xj, X 2 , xj, up to the nth 
point x n . In the same manner, the 2-norm of the xi solution 
trajectory is defined by the positive square root of the sum 
of the squares of absolute values of xi, X 2 , x.?, up to the nth 
point x n . Similarly, the 3-norm for the x solution trajectory 
is defined as the sum of the absolute values (cubed) of xi, 
X 2 , X 3 , up to the nth point x„. The 4-norm for the x solution 
trajectory is defined as the sum of the absolute values (to 
fourth power) of xi, X 2 , x.j, up to the nth point x„. The 5- 
norm for the x solution trajectory is defined as the sum of 
the absolute values (to fifth power) of xi, x 2 , Xj, up to the 
nth point x„. The 6 -norm for the x solution trajectory is 
defined as the sum of the absolute values (to sixth power) of 
Xi, X 2 , xj, up to the nth point x„. The infinity norm for x m is 
defined by the maximum value of the set of the absolute 
values of xi, x 2 , Xi, up to the nth point x„. The same 
procedure holds for the calculation of 1-norm, 2-norm, 3- 
norm, 4-norm, 5-norm, 6 -norm and infinity norm for y and 
z solution trajectories. 

Step V: In the execution program, calculate the difference 
between the solution trajectories of the control model 
equation and the varied model equation by F 1 — x — x m , 
F 2 — y — y m and F 3 — z — z m for the given range of data 
points when j = 1, 2, 3, 4, n such that the difference 
between the solution trajectories of the control model 
equations and the varied model equations would be F 1 = 
Xj x m j, F 2 yj y-mj an d F 3 Zj z m j. 

For the purpose of this analysis, we will also calculate the 
1-norm, 2-norm, 3-norm, 4-norm, 5-norm, 6 -norm and 
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infinity norm of F 1 ,F 2 and F 3 . For example, 1-norm of Fi 
will be the sum of the absolute values of the data 
points (*! - x ml ), (x 2 - x m2 ), (x 3 - x m 3 ),{x A - x m4 ), ..., 
( x n — x mn ) where x 1 and x ml stand for the first data point 
of x solution trajectory and the first data point of the 
modified x m solution trajectory respectively. x 2 andx m2 
stand for the second data point of x solution trajectory and 
the second data point of the modified x m solution trajectory 
respectively and so forth. The 1-norm of F 2 will be the sum 
of the absolute values of the data points (y 1 — y), (y 2 — 

/mzl (y 3 -y m 3)-(y 4 -y m4 ). •••> (y n -ymJ where y x 

and y ml stand for the first data point of y solution trajectory 
and the first data point of the modified y m solution 
trajectory respectively. yandy m2 stand for the second data 
point of y solution trajectory and the second data point of 
the modified y m solution trajectory respectively and so 
forth. Similarly, the 1-norm of F 3 will be the sum of the 
absolute values of the data points [z 1 — z ml ), ( z 2 — z m2 ), 
(z 3 ~ z),(z 4 - z m4 ), ..., (z n -z mn ) where z x and z ml 
stand for the first data point of z solution trajectory and the 
first data point of the modified z m solution trajectory 
respectively. z 2 andz 7n2 stand for the second data point of z 
solution trajectory and the second data point of the modified 
z m solution trajectory respectively and so forth. The 2- 
norm, 3-norm, 4-norm, 5-norm, 6 -norm and infinity norm 
can similarly be calculated for the differences of three 
solution trajectories F 1 ,F 2 and F 3 . 

Step VI: To calculate the cumulative percentage effect of 
variation of a chosen model parameter one-at-a-time when 
other parameters are fixed on each solution trajectory. For x 
solution trajectory, we will calculate the following values: 
(1-norm of F 1 divided by the 1-norm of x) multiplied by 
100; (2-norm of F 1 divided by the 2-norm of x) multiplied 
by 100; (3-norm of F r divided by the 3-norm of x) 
multiplied by 100; (4-norm of F 1 divided by the 4-norm of 

x) multiplied by 100; (5-norm of F 1 divided by the 5-norm 
of x) multiplied by 100; ( 6 -norm of F 1 divided by the 6 - 
norm of x) multiplied by 100 and (infinity-norm of F 1 
divided by the infinity-norm of x) multiplied by 100 . 

To calculate the percentage cumulative effect of variation of 
a model parameter one-at-a-time when other parameters are 
fixed on y solution trajectory, we will calculate the 
following values: (1-norm of F 2 divided by the 1-norm of v) 
multiplied by 100; (2-norm of F 2 divided by the 2-norm of 

y) multiplied by 100; (3-norm of F 2 divided by the 3-norm 
of v) multiplied by 100; (4-norm of F 2 divided by the 4- 
norm of y) multiplied by 100; (5-norm of F 2 divided by the 
5-norm of y) multiplied by 100; ( 6 -norm of F 2 divided by 
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the 6-norm of y) multiplied by 100 and (infinity-norm of F 2 
divided by the infinity-norm of y) multiplied by 100. 

To calculate the percentage cumulative effect of variation of 
a model parameter one-at-a-time when other parameters are 
fixed on z solution trajectory, we will calculate the 
following values: (1-norm of F 3 divided by the 1-norm of z) 
multiplied by 100; (2-norm of F 3 divided by the 2-norm of 
z ) multiplied by 100; (3-norm of F 3 divided by the 3-norm 
of z ) multiplied by 100; (4-norm of F 3 divided by the 4- 
norm of z) multiplied by 100; (5-norm of F 3 divided by the 
5-norm of z.) multiplied by 100; (6-norm of F 3 divided by 
the 6-norm of z) multiplied by 100 and (infinity-norm of F 3 
divided by the infinity-norm of z) multiplied by 100. 

Due to the unstable values of the 1-norm, 2-norm, 3-norm, 
4-norm, 5-norm, 6-norm and infinity-norm specifications, 
we adopt to use a compact value for related percentage 
values of these norms. For example, the cumulative 
percentage value of 1-norm sensitivity in terms of the 
difference in solution trajectories involves the sum of the 1- 
norm values due to F 1 solution trajectory, 1-norm values 
due to F 2 solution trajectory and 1-norm values due to F 3 
solution trajectory. The same procedure can be followed to 
calculate the cumulative percentage values for 3-norm, 3- 
norm, 4-norm, 5-norm, 6-norm and infinity-norm 
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sensitivities in terms of F 1 solution trajectory, F 2 solution 
trajectory and F 3 solution trajectory. 

The cardinal points of sensitivity analysis results 
interpretation are: 

■S the parameter which when varied a little and 
produces the biggest cumulative effect on the 
solution trajectories is called a most sensitive 
parameter. And it is highly uncertain. This by 
implication will attract error in prediction and thus 
require further validation. 

■S the parameter which when varied a little 
andproduces lower or least sensitivity values is 
called a less or least sensitive parameter and 
requires lesser validation. 

Method for Best-fit parameter value selection 

In order to select the best-fit parameter value for each 

model parameter, we recommend that 

■S At 100% variation, the coordinates of the solution 
trajectories have same values, sums, squares and 
square roots = 0. Therefore, the norm values at 
100 % are zero. 

The Local minimum value is selected at a point where the 
smallest norm value occurs before or after the 100% 
variation. 


IV. RESULTS AND DISCUSSION 

On the application of the define method of analysis, we hereby present and discuss the following results: 


r 

1 -norm 

Table. 1: 

2-norm 3-norm 

sensitivity analysis results for model parameter r 

4-norm 5-norm 6-norm co-norm 

0.0030 

143.9520 

78.8284 

82.0093 

66.2747 

56.624 

50.6594 

29.8099 

0.0060 

143.1933 

78.4409 81.6042 

65.9501 

56.348 

50.4130 

29.6754 

0.0090 

142.3997 

78.0321 

81.1766 

65.6070 

56.056 

50.1521 

29.5334 

0.0120 

141.5698 

77.6010 

80.7256 

65.2447 

55.747 

49.8763 

29.3837 

0.0150 

140.7023 

77.1468 

80.2503 

64.8624 

55.421 

49.5850 

29.2259 

0.0180 

139.7958 

76.6685 

79.7499 

64.4593 

55.077 

49.2776 

29.0605 

0.0210 

138.8491 

76.1654 

79.2236 

64.0350 

54.715 

48.9538 

28.8879 

0.0240 

137.8610 

75.6368 

78.6706 

63.5887 

54.334 

48.6131 

28.7078 

0.0270 

136.8304 

75.0818 

78.0903 

63.1201 

53.934 

48.2553 

28.5200 

0.0300 

135.7565 

74.4999 

77.4821 

62.6287 

53.514 

47.8800 

28.3243 


What do we learn from Table 1? 


We have observed from table 1 that as the value of the 
model parameter r increases monotonically from 0.003 
(approx) to 0.03 (approx), the 1-norm error data decreases 
monotonically from the value of 143.95 (approx.) to 135.76 
(approx.), the 2-norm error data decreases monotonically 
from the value of 78.83 (approx.) to 75.50 (approx.), the 3- 
norm error data decreases monotonically from the value of 
82.01 (approx.) to 77.48 (approx.), the 4-norm error data 


decreases monotonically from the value of 66.27 (approx.) 
to 62.63 (approx.), the 5-norm error data decreases 
monotonically from the value of 56.62 (approx.) to 53.51 
(approx.), the 6-norm error data decreases monotonically 
from the value of 50.66 (approx.) to 47.88(approx.) and the 
infinity-norm error data decreases monotonically from the 
value of 29.81 (approx.) to 28.32 (approx.). 
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On the basis of this present analysis we have observed that r 
= 0.003 is associated with relatively higher uncertainty 
when compared to the value of r = 0.03 irrespective of the 
type of p-norm we have used to calculate the sensitivity 
analysis. Despite the observed uncertainty analysis results 
due to a one percent to ten percent variation of the intrinsic 
growth rate of the fish population, it is clear that the 
statistical range of the p-norm sensitivity values are listed as 
follows: the range of 1-norm statistical range is 8.19, the 
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range of 2-norm statistical range is 3.33,the range of 3-norm 
statistical range is 4.53,the range of 4-norm statistical range 
is 3.64,the range of 5-norm statistical range is 3.11,the 
range of 6-norm statistical range is 2.78 .the range of 
infinity-norm statistical range is 1.49. 

These bigger sensitivity values which indicate high 
uncertainty of the intrinsic growth rate of the fish 
population can be further minimized. 


Table.2: local minimum resultsfor parameter r 

r 1-norm 2-norm 3-norm 4-norm 5-norm 6-norm oo-norm 


0.2880 

4.3024 


2.2917 

2.4375 2.0011 


1.739 

1.5812 1.0897 

0.2910 

3.2128 

1.7106 

1.8194 

1.4936 

1.298 

1.1803 

0.8137 

0.2940 

2.1326 

1.1350 

1.2072 

0.9910 

0.861 

0.7831 

0.5401 

0.2970 

1.0617 

0.5648 

0.6007 

0.4931 

0.429 

0.3897 

0.2688 

0.3000 

0.0000 

0.0000 

0.0000 

0.0000 

0.000 

0.0000 

0.0000 

0.3030 

1.0527 

0.5595 

0.5951 

0.4885 

0.425 

0.3860 

0.2663 

0.3060 

2.0964 

1.1138 

1.1847 

0.9724 

0.845 

0.7685 

0.5304 

0.3090 

3.1314 

1.6630 

1.7688 

1.4517 

1.262 

1.1473 

0.7926 

0.3120 

4.1577 

2.2071 

2.3475 

1.9265 

1.674 

1.5225 

1.0527 

0.3150 

5.1755 

2.7463 

2.9209 

2.3969 

2.083 

1.8943 

1.3106 


What do we learn from Table 2? 

Table 2 shows a result of the cumulative effect of 96 to 105 
percent variation of model parameter r. At lOOpercent 
variation, there are no changes in the original and the varied 
solution trajectories, hence the zero values. 

Observe that as the model parameter value increase 
monotonically from a low value of 0.2880 to 0.0000 
corresponding to 100 percent variation and then increases 
monotonically to 0.3150 approximately in column I, the 1- 
norm sensitivity value decreased monotonically from a 
value of 4.30 to 0.0000 corresponding 100 percent variation 
and then increases monotonically to 5.18 approximately in 
column II, the 2- norm sensitivity value decreased 
monotonically from a value of 2.29 to 0.0000 corresponding 
lOOpercent variation and then increases monotonically to 
2.75 approximately in column III, the 3- norm sensitivity 
value decreased monotonically from a value of 2.44 to 
0.0000 corresponding lOOpercent variation and then 
increases monotonically to 2.92 approximately in column 
IV, the 4- norm sensitivity value decreased monotonically 
from a value of 2.00 to 0.0000 corresponding lOOpercent 
variation and then increases monotonically to 2.40 
approximately in column V, the 5- norm sensitivity value 
decreased monotonically from a value of 1.74 to 0.0000 
corresponding lOOpercent variation and then increases 
monotonically to 2.08 approximately in column VI, the 6- 


norm sensitivity value decreased monotonically from a 
value of 1.58 to 0.0000 corresponding lOOpercent variation 
and then increases monotonically to 1.90 approximately in 
column VII, the infinity- norm sensitivity value decreased 
monotonically from a value of 1.09 to 0.0000 corresponding 
lOOpercent variation and then increases monotonically to 
1.31 approximately in column VIII. 

The local minimum is selected at the greatest lower bound 
of the data points for each p-norm results. For instance, the 
local minimum value for model parameter r is 0.303 where 
1-norm local minimum value = 1.0527, 2-norm local 
minimum value = 0.5595, 3-norm local minimum value = 
0.0.5951, 4-norm local minimum value = 0.4885, 5-norm 
local minimum value = 0.4250, 6-norm local minimum 
value = 0.3860 and infinity-norm local minimum value = 
0.2663. 

V. CONCLUSION 

This present study has been able to reduce the uncertainty in 
the intrinsic growth rate of a fish population due to its 
variations using a combination of a numerical scheme and 
the mathematical p-norms. The present results compliment 
the earlier contribution of Ekaka-a et al (2012) that only 
considered the sensitivity analysis in the context of a shorter 
experimental time. This present proposed numerical scheme 
can be extended to study the sensitivity analysis of other 
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model parameter values which we did not consider in this 
study in our future investigation. 
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